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The conditions under which quantum-classical Liouville dynamics may be reduced to a master 
equation are investigated. Systems that can be partitioned into a quantum-classical subsystem 
interacting with a classical bath are considered. Starting with an exact non-Markovian equation for 
the diagonal elements of the density matrix, an evolution equation for the subsystem density matrix 
is derived. One contribution to this equation contains the bath average of a memory kernel that 
accounts for all coherences in the system. It is shown to be a rapidly decaying function, motivating a 
Markovian approximation on this term in the evolution equation. The resulting subsystem density 
matrix equation is still non-Markovian due to the fact that bath degrees of freedom have been 
projected out of the dynamics. Provided the computation of non-equilibrium average values or 
correlation functions is considered, the non-Markovian character of this equation can be removed by 
lifting the equation into the full phase space of the system. This leads to a trajectory description of 
the dynamics where each fictitious trajectory accounts for decoherence due to the bath degrees of 
freedom. The results are illustrated by computations of the rate constant of a model nonadiabatic 
chemical reaction. 



I. INTRODUCTION 

When investigating quantum relaxation processes in 
the condensed phase, one often partitions the full quan- 
tum system into a subsystem whose dynamics is of in- 
terest and an environment or bath with which the sub- 
system interacts. There is a large literature dealing with 
such open quantum systemsi^ A number of different 
equations of motion for the density matrix of the sub- 
system have been derived, including the Lindblad^ and 
Redfield^ equations and a variety of generalized quantum 
master equations 1 1 1 1 Effects due to the environ- 

ment typically enter these equations through coupling 
terms involving parameters that characterize the bath 
relaxation processes. Such equations have been used to 
investigate aspects of decoherence in the quantum sub- 
system arising from interactions with bath degrees of 
freedom. 

Sometimes it is convenient to suppose that the dynam- 
ics of certain degrees of freedom is described by quantum 
mechanics while other degrees of freedom may be treated 
by classical mechanics to a good approximation. This is 
the case if one considers systems involving light particles 
interacting with a bath of heavy particles. Proton and 
electron transfer processes in the condensed phase and 
in biomolecules fall into this category, as do many vibra- 
tional relaxation processes. One is then led to study the 
dynamics of quantum-classical systems where the entire 
system is partitioned into quantum and classical subsys- 
tems. 12 Equations of motion for the quantum subsystem 
density matrix, where the classical bath is modeled as a 
dissipative environment have been derived i 13 i 14 ' 15 Such 
descriptions are useful for many applications; however, 
there are situations where the quantum subsystem evo- 
lution depends explicitly on the details of the bath dy- 
namics. This is important since specific features of bath 
motions can influence quantum rate processes^ To de- 



scribe such specific bath dynamical effects one must use 
the full quantum-classical equation of motion. 

Another partition of the system is required when the 
quantum subsystem is directly coupled to a small subset 
of the environmental degrees of freedom. For example, in 
proton transfer within a bio-molecule, the quantum sub- 
system may be taken to be the proton, which interacts 
directly with a specific set of functional groups of a larger 
molecule immersed in a solvent. In this case, we may 
define a quantum-classical subsystem comprising both 
quantum degrees of freedom (the proton) and a subset of 
the classical variables that directly couple to these quan- 
tum degrees of freedom (the specified functional groups). 
The remaining classical variables constitute the bath. In 
such a partition, the bath may be treated either explicitly 
or as a dissipative environment. Equations of motion for 
the density matrix of a quantum-classical subsystem in- 
teracting with a dissipative bath have been derived^ In 
this article we consider quantum-classical systems of this 
type but instead retain the details of the bath dynamics 
and study the conditions under which the dynamics can 
be reduced to a master equation. 

The reduction of quantum-classical dynamics to a sub- 
system master equation hinges on the decoherence in the 
subsystem induced by interactions with the bath degrees 
of freedom^ Consequently, we focus on how decoher- 
ence is described in quantum-classical systems and the 
conditions under which it is strong enough to eliminate 
the off-diagonal elements of the density matrix on short 
time scales. The result of the analysis is a non-Markovian 
generalized master equation for the density matrix of the 
quantum-classical subsystem. When this equation is used 
to compute non-equilibrium averages or correlation func- 
tion expressions involving subsystem properties, we show 
that the subsystem dynamics may be lifted to the full 
phase space, including the bath degrees of freedom, to 
recover a Markovian master equation for the quantum 
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and all classical degrees of freedom. Solutions of this 
equation may be obtained from an ensemble of surface- 
hopping trajectories, each member of which incorporates 
the effects of quantum decoherence. 

The outline of the paper is as follows: The starting 
point of our analysis is the quantum-classical Liouville 
equation for the entire system ^^1^2^^^^.^^^^ 
In Sec. [n] we show how this equation can be cast into the 
form of a generalized master equation for the diagonal 
elements of the density matrix. This equation involves 
a memory kernel operator that contains all information 
on quantum coherence in the system. We discuss the ex- 
plicit form of the memory kernel operator that governs 
the evolution in off-diagonal space and present its func- 
tional form. The analysis in Sec. IIIII and Appendix [B] 
shows that an average over an ensemble of trajectories 
evolving on coherently coupled surfaces with different 
bath initial conditions can be used to reduce the full- 
system generalized master equation to a non-Markovian 
subsystem generalized master equation. We show that 
this equation can then be lifted back to the full phase 
space to obtain a Markovian master equation. We apply 
this formalism in Sec. HVI and calculate the nonadiabatic 
rate constant for a model system. Numerical results us- 
ing full quantum-classical Liouville dynamics and mas- 
ter equation dynamics are compared. In the concluding 
section we comment on the relationship of our master 
equation dynamics to other surface-hopping methods. 



II. GENERALIZED MASTER EQUATION 

Starting from the quantum-classical Liouville equa- 
tion, it is not difficult to derive a generalized master equa- 
tion for the diagonal elements of the density matrix. As 
described in the Introduction, the position and momen- 
tum operators, {q,p), of the quantum degrees of freedom 
are assumed to be coupled directly to a set of classical 
phase space variables, X = (R ,P ). Together these 
make up the quantum-classical subsystem. The classical 
Xq variables are, in turn, directly coupled to the remain- 
der of the classical phase space variables, Xb = (Rb, Pb), 
that constitute the bath. (Our formulation must be mod- 
ified If the quantum degrees of freedom couple directly to 
all other variables in the system.) The total Hamiltonian 
of the system is, 
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The potential energy operator, V(q, Ro, Rb), includes the 
contributions from the quantum-classical subsystem, the 
bath, and the interaction between the two. It is often 
convenient to represent the dynamics in the adiabatic 
basis given by, h(R 0l R b )\a; R ) = E a (R , R b )\a; R ), 
where h(Ro, Rb) is the quantum Hamiltonian for a fixed 



configuration of the classical particles. The adiabatic 
eigenf unctions depend only on the coordinates Rq since 
the dependence on the bath coordinates enters through 
the potential as an additive constant. In this basis, the 
equation of motion for the full density matrix is^, 



dt 



Pw (X,t) 
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where py^* (X, t) is a density matrix element which de- 
pends on the full set of phase space coordinates X — 
(R,P) = (Xo,Xb). The quantum-classical Liouville 
super-operator is given by^ 

— i£ a a',00' = —ii^aa' + L aa ')S a jsd a ' pi + Jaa',00' , (3) 

where uj aa > = AE aa ' /H, with AE aa < = E a — E a i , and 
the classical Liouville operator, iL aa i is given by 
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iL aa , = — . — + -[F c 
M dR 2 



F° 



The Hellmann-Feynman forced for state a is 
F a = {a- 1 R \dV(q : Ro,Rb)/dR\a;Ro), and the op- 
erator Saa', 00' i defined in the next section, accounts 
for quantum transitions and corresponding momentum 
changes in the environment. 

We shall often simplify the notation in what fol- 
lows. Above, pw(X,t) refers to the partially Wigner 
transformed density matrix whose matrix elements are 
p'ff (X,t). Since we use partially Wigner transformed 
variables throughout this article, we shall drop the sub- 
script W. 

The quantum-classical Liouville evolution operator can 
be partitioned into diagonal, off-diagonal and coupling 
components by defining the super-operators: £ d , £ do , 
C°' d , and £°, where the d, and o superscripts denote 
diagonal and off-diagonal, respectively. We define the di- 
agonal part of the density matrix p d (X,t), with matrix 
elements, p aa (X, t)S aa > = p°t(X,t). Similarly, the off- 
diagonal part of the density matrix, p a (X,t), has matrix 
elements p aa ' (X,t){l - 8 aa >) = p% a ' (X,t). Using these 
definitions, the quantum-classical Liouville equation may 
be expressed formally as the following set of coupled dif- 
ferential equations, 



_5 

dt 



Pd {X,t) = -i£ d Pd (X,t)~i£ d >°p (X,t) (5) 
Po (X,t) = -i£°p (X,t)-i£ od p d (X,t). (6) 



By substituting the formal solution of Eq. © into 
Eq. (O, we obtain the evolution equation for p d (X,t), 



-i£ d '°e' lCOt Po (X, 0) - i£ d p d (X, t) (7) 



dt'iC d < e- iC0{ - t - t '}i£ ' d Pd {X,t') 



For the remainder of this analysis we will assume that 
p o (X,0) — 0, and thus the first term vanishes. This 
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amounts to initially preparing the system in a pure state 
or incoherent mixture of states. Although equation ([7]) 
is general and may be used to study systems that are 
initially prepared in coherent states, we shall not con- 
sider such situations here. Using Eq. ([3]), the explicit 
definitions of the matrix elements of the Liouville super- 
operators are, 



iC d 

lL 'aa',00' 



iC 
iC 



act' ,00' 
o . d 



(8) 



= iC aa ',00'V- — 5 aa >)(l — (5/3/3'); 



act', 00' 



iC° 



,00' 



where iL a = iL aa and J^ fj , = J aa ,00'> for (3 ^ 0', with 
a similar definition for J°' d . Using these definitions and 
the initial condition discussed above, we obtain the gen- 
eralized master equation for the evolution of the diagonal 
elements of the density matrix. 
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p a d (X,t) = -iL a fi(X,t) 



+ I dt'J2M a 0(t')p^X,t-t'),(9) 
Jo 

where the memory kernel operator A4 a p(t) is given by, 
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(10) 

and acts on all of the classical degrees of freedom that ap- 
pear in functions to its right. Next, we analyze the form 
of the memory kernel operator (|10p in order to cast it 
into a form that is suitable for the derivation of a master 
equation. 



Memory kernel 

The explicit form of the J operator was derived pre- 
viously 33 and is given by 



Jaa',00' 



C a 0d a '0' + C%/ a /6, 



a '0'°a0, 



(11) 



where 



C a p = -D afi {X ) (l + ~S a/3 ■ JL} , (12) 

the nonadiabatic coupling matrix element is given by 
d a = (a; i?o|V r \(3; Ro) , S a = AE a pd a 0l 'D a p(Xo), 
and D a p{X ) = (P /M ) ■ d a p. We have shown in earlier 
work that the action of the C operator on phase space 
functions may be computed using the momentum-jump 
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approximation 

C a 0(Xo) 



-D a 0(Xo)j a 0(X o ) 



(13) 



where the momentum shift operator, j a 0(Xo), is a trans- 
lation operator in momentum space, 



j a0 f(P O ) = e ^M o9/ 9(P .d a ^ f{Po) 
= f(P O + AP Oa0 ), 

APa aj3 = d a (sgn(P • d a p) 



(14) 
(15) 



x V (Flo ■ d a0 ) 2 + AE a/3 M - (P ■ d aP ) 



Since momentum shifts occur in conjunction with quan- 
tum transitions, they depend on the quantum states in- 
volved in the transition. Consequently, we use the follow- 
ing notation: X 0af3 = (R Q ,_P Oa 0) = (Ro, Po + AP Oa0 ), so 
that, j a p(X )f(X ) = f(X 0af3 ). It is worth noting here 
that the momentum shift operators do not act on the 
full classical environment. They only act on the classical 
variables Xq that directly couple to the quantum degrees 
of freedom. We also observe that the argument of the 
square root in Eq. (|15[) must be positive. This condition 
prevents quantum transitions when there is insufficient 
energy in the classical degrees of freedom to effect the 
transition. 

The time evolution in Eq. (|10p , is given by the propaga- 
tor, e — The quantum-classical Liouville operator 
in this expression acts on the entire phase space X and 
accounts for the following processes: classical evolution 
of the bath coordinates X b and evolution of the classical 
subsystem coordinates Xq on the mean potential surface 
(E^ + E' fl )/2 with an associated phase factor. This evolu- 
tion is interspersed with quantum transitions taking the 
subsystem to other coherently coupled states where evo- 
lution is again on mean surfaces with associated phase 
factors. In the course of this evolution the system never 
returns to a diagonal state involving evolution on a sin- 
gle adiabatic surface. For this reason we refer to such 
evolution as being in "off-diagonal space" . 

If one considers the explicit action of the propagators 
appearing in the memory kernel, one can show how this 
operator acts on an arbitrary function of the phase space 
coordinates X . The details are given in Appendix [A] 
Using the results obtained there, one can show that the 
generalized master equation can be written as, 
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at 



fi(X,t) = -iL a fi(X,t) 



(16) 



dt' ( £ M${X, (X^ t , , X bit , , t -t') 



■ Mgg (X, t')pd (X^ v t , , X bit > , t - t 1 ) 



In this expression the superscripts on the coordinates in- 
dicate the action of a second momentum shift operator, 
(j va (X au)f(X au,X b ) = f(X^,X b j). This result pro- 
vides us with a definition of the memory function, 



Mf g (X,t) = 2Rc 



L u0 



W a 0(t', 0)J D aP (X Q )D a p(X w ,), 

(17) 
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where the subscripts and superscripts on the memory 
function label the indices on the first and second D func- 
tion respectively. The phase factor W a p is defined as 

W a0 (h,t 2 ) = e~ lf t? ^"^(-fW.r) . (i 8 ) 

Now the actions of all classical propagators and momen- 
tum jumps have been accounted for explicitly in Eq. (|16[) 
so that the memory kernel is a function of the phase space 
variables. 

Given the assumption about the initial condition on 
the density matrix (p o (X,0) — 0), for a two-level sys- 
tem the generalized master equation (fT6|) is fully equiv- 
alent to the quantum-classical Liouville equation from 
which it was derived. This result is also applicable to 
multi-level systems in a weak coupling limit where quan- 
tum transitions among different coherently coupled states 
are neglected in the the off-diagonal propagator. This 
amounts to neglecting terms higher than quadratic order 
in the nonadiabatic coupling strength, d a /3, in the evolu- 
tion operators. The time evolution described by Eq. (fl6|) 
consists of classical evolution along single adiabatic sur- 
faces and two memory terms. The memory terms account 
for transitions to the mean surface, evolution along this 
surface, and transitions to a new adiabatic surface with 
rate M™p, or transitions back to the original surface with 
rate M^". Thus, the dynamics of the generalized mas- 
ter equation derived here is separated into diagonal and 
off-diagonal components providing a framework within 
which to investigate decoherence in the quantum-classical 
subsystem induced by the bath. 



III. MASTER EQUATION 

We next consider the conditions under which the gen- 
eralized master equation (|16p may be reduced to a sim- 
ple master equation without memory. This reduction 
hinges on the ability to consider the memory kernel as a 
rapidly decaying function so that a Markovian approxi- 
mation can be made. From its form in Eq. (|17[) one can 
see that M(X, t), which contains all information on quan- 
tum coherence, is an oscillatory function. As a result, a 
Markovian approximation to the memory kernel cannot 
be made directly on the full phase space equation since 
there is no obvious mechanism for the decay of the mem- 
ory function. It is the decoherence by the environment 
that provides such a mechanism. 

In this analysis we exploit the fact that decoherence 
has its origin in interactions with the bath degrees of 
freedom. We have already observed that we are inter- 
ested in dynamical properties of the quantum-classical 
subsystem. For instance, nonequilibrium average values 



of interest have the form, 




= W dX A^(X )pf(X ,t) , (19) 

a/3 

where A /3a (Xo) are the matrix elements of a property of 
the subsystem and p^(X ,t) = J dX b p aP {X,t) is the 
subsystem density matrix. If the operator A /3a (Xo) is 
diagonal, then only the diagonal elements of the subsys- 
tem density matrix are needed to compute its average 
value. Alternatively, if decoherence quickly destroys the 
off-diagonal subsystem density matrix elements, then, af- 
ter a short transient, only the diagonal elements will be 
needed to compute the expectation value. Later, we shall 
show that similar considerations can be used to evaluate 
correlation function expressions for transport properties 
of the subsystem. 

To compute such average quantities, we see that we 
need the subsystem density matrix elements. Starting 
with the generalized master equation in full phase space 
(Eq. (|16|) ). we can introduce a bath projection opera- 
tor^ 

V- = p c b {X b ;R ) JdX b - , (20) 

whose action on the density matrix yields the subsystem 
density matrix. Here pl(X b ; Rq) is the bath density con- 
ditional on the subsystem configuration space variables 
Rq. The bath density function is in general quantum 
mechanical but in some applications it may be replaced 
by its high temperature classical limit. In the projection 
operator formalism it is not necessary to distinguish be- 
tween these two cases. The complement of V is Q. In 
Appendix [B] using standard projection operator meth- 
ods^ we argue that bath averaged correlations involv- 
ing the fluctuations of the memory kernel from its bath 
average may be neglected. In this case, the subsystem 
generalized master equation is, 

^(X ,t) = (21) 
-(tL a e- iQL ^Q P 2(X,0)) b - (iL a ) b p™(X ,t) 

- f dt'{iL a e- iQL ^iQL a ) h p a s {X ,t-t') 
Jo 

Jo v p 
+ J2(M^(X,t')) b p a d (X^ t „t- o) . 

V 

In this expression the memory function appears in the 
form of an average over a bath distribution function con- 
ditional on the subsystem configuration space variables, 
(M^(X,t')) b = f dX b M^{X,t')p c b (X b ;Ro). Since the 
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memory function M^(X,t') involves evolution on the 
mean of the a and j3 adiabatic surfaces, the average over 
bath initial conditions will result in an ensemble of such 
trajectories, each carrying an associated phase factor. As 
a result, the bath ensemble average (M"^(X,t')) b will 
decay on a time scale characterized by the decoherence 
time, Tdecoh- If the decoherence time is short compared 
to the decay of the populations we can make a Markovian 
approximation, 

(M^(X,t')) b « 2(j™ dt'{M«l{X,t'))^5{t') 

= 2m af} (X )5(t') . (22) 

Applying this Markovian approximation to Eq. (|21| . we 
obtain, 

§- t P«(Xo,t)= (23) 
~ J dX b iL a e~ tQL ^Q P a d {X,Q) - (iL a ) b p«(X ,t) 

- f dt'iiL^^'iQLMiX^t-t') 

+ ^2m a p(X )j a ^ l3 pP(X ,t) - m aa (X )Ps{X ,t) , 

P 

where 

poo 

m aa (Xo) = ~ J2 / dt'(M™(X,t')) b . (24) 

The use of the Markovian approximation results in the 
instantaneous action of two momentum shift operators 
on the density. Consequently, the penultimate term in 
Eq. (|23p was rewritten to incorporate a single momentum 
shift operator whose action is j a -,/3f(X ) — f(Xg^) re- 
sulting from a transition from one single adiabatic surface 
to another single adiabatic surface: 

ja^f(Po) = j a p(X )j a p(X )f(Po) 

= /(P + AP ^), (25) 

& p Cp = 4/3(sgn(P • 4/3) (26) 

x yj (P ■ d aP ) 2 + 2AP Q/3 M - (P • 4/3)) ■ 

This momentum shift differs from that defined earlier by 
a factor of 2 in front of the energy difference since this 
operator captures the action of two jumps. In the last 
term in Eq. (j23|) the net effect of two momentum shift 
operators with reversed indices acting simultaneously is 
jva(Xoav)jav(X )f(X) = f(X). Since these momentum 
shift operators are inverses of each other, there is no net 
shift. 

As discussed in Appendix [XI the transition rate 
M™p(X, t) captures the effect of two momentum shift op- 
erators. Each of these operators imposes a condition on 



the subsystem kinetic energy that ensures transitions to 
and from the mean surface are allowed. Consequently, 
the transition rate m a /3(Xo) inherits these conditions. 
For example, if a < [3 then m a p(Xo) is non-zero only 
if (P • d a p) 2 /2M > AE fja . Conversely, if a > /3 there 
is no such restriction. In contrast, the transition rate 
m aa (Xo) is non-zero only if (P • d av ) 2 /2M > 
This condition arises from the fact that this contribution 
has its origin from transitions to the mean surface and 
then back to the original surface. 

Lift to full phase space 

Equation (f2"3"|) is still rather difficult to solve since it 
contains a convolution involving bath projected dynam- 
ics. Often, the non-Markovian character of an equation 
can be removed by expanding the space upon which the 
equation is defined. In the analysis above, the non- 
Markovian character arose by projecting out the bath 
variables to obtain a description in the subsystem phase 
space. Consequently, by lifting this equation back into 
the full phase space we can recover the Markovian nature 
of the dynamics. In the full phase space the equation of 
motion is given by 

^fi(X,t) = -iL a fi(X,t) (27) 

+ ^2m aP (X )j a ^pp l3 d (X,t) -m aa (X )p d t (X,t) . 



It is easily verified that applying the projection operator 
algebra to this equation, using the projection operator 
defined in Eq. (f20|) , one obtains the subsystem evolu- 
tion equation (I23|) . Thus, when computing average val- 
ues like those in Eq. (|19|) . or their correlation function 
analogs discussed below, the master equation lifted to 
full phase space yields results identical to those of the 
non-Markovian equation (|23[) . 

Through this analysis we have succeeded in finding a 
master equation description of the dynamics in the full 
phase space which incorporates the effects of decoher- 
ence. The first term in Eq. (|2T|) yields dynamics on sin- 
gle adiabatic surfaces. The other terms correspond to 
contributions to the evolution due to nonadiabatic tran- 
sitions between adiabatic states. The nonadiabatic tran- 
sition rates in these terms incorporate the effects of de- 
coherence. Given the nature of the dynamics generated 
by this master equation, there is a close connection to 
many currently-used surface-hopping schemes which will 
be discussed below. 

It is instructive to compare master equation and 
quantum-classical Liouville dynamics. The master equa- 
tion (|2T|) . like the full quantum-classical Liouville equa- 
tion (HI), can be simulated by following an ensemble of 
surface-hopping trajectories. The trajectories that enter 
in each description are shown in Fig. [TJ We see that 
in full quantum-classical Liouville dynamics the system 
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makes transitions between single adiabatic surfaces via 
coherently coupled off-diagonal states. Coherence is cre- 
ated when such an off-diagonal state is entered and is 
destroyed when it is left. The average over the ensemble 
accounts for net destruction of coherence in the system 
as it evolves. In contrast, the master equation evolves 
the classical degrees of freedom exclusively on single adi- 
abatic surfaces with instantaneous hops between them. 
Transitions from a diagonal state to a coherently coupled 
state and then back to the diagonal state, which play an 
important role in quantum-classical Liouvillc dynamics, 
are accounted for explicitly in master equation dynamics 
by m aa (Xo). Each single (fictitious) trajectory accounts 
for an ensemble of trajectories that correspond to differ- 
ent bath initial conditions. In this connection the evo- 
lution in off-diagonal space is crucial: for a given initial 
subsystem coordinate, the choice of different bath coor- 
dinates will result in different trajectories on the mean 
surface. Thus, it is the average over this collection of 
classical evolution segments that results in decoherence. 
Consequently, this master equation in full phase space 
provides a description in terms of fictitious trajectories, 
each of which accounts for decoherence. When the ap- 
proximations that lead to the master equation are valid, 
this provides a useful simulation tool since no oscillatory 
phase factors appear in the trajectory evolution. 




b) 



X. 



X. 



FIG. 1: Trajectories that enter the solution of the quantum- 
classical Liouville and master equations, (a) In quantum- 



classical Liouville dynamics, the system makes transitions be- 
tween single adiabatic surfaces and coherently coupled states 
involving evolution on mean surfaces. |(b)| In master equa- 
tion evolution, the dynamics is restricted to single adiabatic 
surfaces. All off-diagonal evolution is accounted for by the 
memory kernel. 



IV. APPLICATION TO REACTION RATES 



In this section we apply the formalism developed above 
to calculate the rate constants of a reaction A ^ B. 
For this reaction, the quantum-classical forward rate con- 
stant was derived earlier and is given by 40 



kAB (t) 



dXRe [N% a ' {X, t)W2' a (X, ^)] , 



X 



where N^ a (X, t) is the time evolved matrix element 
of the number operator for the product state B. At 
t = this operator is diagonal in the adiabatic basis 
and its value, in this context, depends only on the sub- 
system coordinates Rq. The spectral density function, 
W2 a (X,ihf]/2), accounts for the quantum equilibrium 
structure of the entire system i 40 i 41 The spectral density 
can be approximated by the form W% a (X, ih(3/2) as 
W2' a (X ,ih/3/2)pl(X b ; R ), such that it is factorized 
into subsystem and bath components. 42 Performing the 
integration over the bath variables, the rate constant ex- 



pression may be written as, 



kAs(t) = 



1 



dX dX b N% a (X,t)pftX h ;Ro) 

xW2 a (X ,^) (29) 
./.Y„Ro (N° a '(X,t)) b W%' a (X ,^) 



We see that the calculation of the rate coefficient entails 
knowledge of the bath average of the time-evolved species 
variable and sampling from the subsystem spectral den- 
sity function. The time evolution of this species variable 
may be calculated using mixed quantum-classical dynam- 
ics^ In general, the subsystem spectral density contains 
both diagonal and off-diagonal components; therefore, 
both diagonal and off-diagonal components of the species 
operator contribute to the computation of the rate coef- 
ficient. Previous work has shown that the off-diagonal 
contributions to the rate constant are negligible^ allow- 
ing one to consider only diagonal contributions. 

The computation of the time evolution of the bath av- 
eraged species variable is completely analogous to the 
calculation of the subsystem density matrix leading to 
Eq. (|23[k however, now the analysis must be carried out 
starting with the quantum-classical Heisenberg equation 
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of motion, 33 



Strong Coupling 



Weak Coupling 



-A« (X,t) 



iC aa , M ,AW'(X,t). (30) 



The rate coefficient can be computed from the expression 
in the first line of Eq. (f29|) using the lifted form of the evo- 
lution equation for the diagonal elements of a dynamical 
variable, 



dt 



A aa (X,t) = iL a (X )A aa (X,t) 



(31) 





where the memory function, m\ is the adjoint of m de- 
fined previously in Eq. |2"2"1) . The effects of decoherence 
that lead to this expression restrict the evolution of the 
observable to its diagonal components. Therefore, one 
only needs to consider the diagonal terms of the subsys- 
tem spectral density in the calculation of the rate coeffi- 
cient. 



Model system 

As an application of this formalism, we consider a sim- 
ple model for a quantum rate process that has been stud- 
ied earlier using mixed quantum-classical dynamics. 44 
The investigation of this model allows us to assess the 
validity of the Markovian approximation, Eq. (|22p . and 
the utility of the master equation for calculating the rate 
coefficient 

The model is a two-level system to which we couple 
v oscillators. The subsystem consists of the two-level 
quantum system bilinearly coupled to a non-linear oscil- 
lator with phase space coordinates (Ro, Pq) governed by 
a symmetric quartic potential, V q (Ro) — aR^/A — bR^/2. 
The bath consists of v — 1 = 300 harmonic oscillators 
whose frequencies u>j, are distributed with Ohmic spec- 
tral density that depends on the Kondo parameter^ 
The bath is bilinearly coupled to the subsystem oscillator 
such that the quantum system does not directly interact 
with the bath; it only feels its effects through the cou- 
pling to the quartic oscillator. As discussed elsewhere, 
it has been argued that this model captures many of the 
essential features of condensed phase proton transfer pro- 
cesses 

Using a diabatic representation, the Hamiltonian for 
this system is, 



H 



+ 



V q (R ) + h l0 R -Ml 

-m v q {R ) - h 10 Ro 



P 2 
2M n 



E 



It 

i 3 



R i 




(32) 




FIG. 2: Plots of free energy vs Ro for strong and weak cou- 
pling cases. The parameters are: 70 = 10.56 for strong cou- 
pling and 70 = 2.64 for weak coupling between the two-level 
quantum system and the quartic oscillator. The other pa- 
rameters are the same for both cases: £1 = 0.51, f3 = 0.5, 
£k = 2, A = 0.5 and B = 1. These parameters were chosen 
to give a well-defined rate process with a significant num- 
ber of re-crossing events. The small energy gap ensures that 
the majority of trajectories satisfy the energetic requirements 
for nonadiabatic transitions given by (|15p and (|26p . and the 
parameter j3 was chosen to be small enough to satisfy the 
high temperature approximation. All other parameters in the 
Ohmic spectral density are the same as those used in earlier 
studies, 44 and the results are presented in the same dimen- 
sionless units as those used in previous studies.—. 



where 2f2 is the adiabatic energy gap. The adiabatic f ree 
energy surfaces, W a (Ro) = V q (Ro)^fh^Q 2 + (70-R0) 2 are 
sketched in Fig. [2] In this figure we also show the mean 
free energy surface, W 12 (R ) = {W l + W 2 )/2 = V q (R ), 
which plays an essential role in the calculation of the 
memory function. 

The simulations of quantum-classical Liouville dynam- 
ics were carried out using the sequential short-time prop- 
agation algorithm 4 ^ in conjunction with the momentum- 
jump approximatio n 20 ! 37 and a bound on the observ- 
able. 20 The initial positions and momenta of the quar- 
tic oscillator and bath were sampled from the classical 
canonical density function. The details of these methods 
can be found elsewherei 20 i 44 ' 46 The simulations of the 
master equation consist of two parts which we describe 
below. First we compute m a /3{Xo) in an independent cal- 
culation involving evolution on the mean surface. Then 
we use this result in the sequential short-time propaga- 
tion algorithm restricted to single adiabatic surfaces. 



Calculation of m a p(Xo) 



In order to investigate the validity of the Markovian 



The solution of the eigenvalue problem for this Hamilto- 
nian yields the adiabatic eigenstates, a;i?o), and eigen- 
values E a (R) = V q (R Q ) + Vb(Rb\ Ro) T ^V^^+T^^o) 2 ' 01 time. From Eq. (f2"§)) . this average is weighted by 



approximation we calculate (M"a(X,t))b, as a function 
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2 0.5 1 1.5 2 2.5 3 
t 



FIG. 3: Plot of the bath averaged memory function 
(Mil (X, i))b versus time for 70 = 2.64, solid line: Ro — 
-0.55, P = 3.2, dotted line: R = 0.4, P = 2.4, dashed line: 
Ro = -0.25, P = -3.4, dot-dash line: R = 0.6, P = -2.2. 
Here we see that for a range of choices of Xo, this function 
decays quickly. 



pI(Xi,;Rq), the Wigner representation of the quantum 
bath distribution conditional on the subsystem coordi- 
nate. In general the determination of the quantum dis- 
tribution function is a difficult problem; however, it is 
known for a harmonic bath^i and may be used to ac- 
count for quantum bath effects. In the Conclusions we 
comment briefly on quantum bath effects in our formal- 
ism. In our calculations we use the high temperature 
limit where the classical canonical equilibrium density, 
conditional on the subsystem configuration, provides a 
good approximation^! 

The quantity (M"^(X,t))b involves the product of the 
initial value of D a j3, the phase factor W a p, and D a p at 
a time-evolved phase point. The latter two quantities 
may be obtained from adiabatic dynamics on the mean 
surface for a given Xq. The bath averaged memory func- 
tion, (M"p(X,t))b may be computed from an average 
over an ensemble of trajectories, each with a fixed ini- 
tial value of Xq and bath coordinates drawn from the 
phase space distribution p%{X\,\ Ro). As discussed above, 
the bath average of this oscillatory function provides a 
mechanism for its decay, characterized by the decoher- 
ence time, Tdecoh- This time will depend on the subsys- 
tem coordinate Xo. In Fig. [3] we plot (M"g(X,t))h as a 
function of time for several subsystem coordinate values 
and show that the bath averaged memory function does 
indeed decay on a rapid time scale. Figure Q] shows how 
the decoherence time, taken as the first zero crossing of 
(M®p(X, t))b, depends on the phase space coordinate Xq. 
In the allowed phase space regions, the decoherence time 
is a relatively weak function of the phase space coordi- 
nates, with the exception of some localized regions where 
it varies strongly. From these results we may compute 
the mean decoherence time and find Tdecoh 

= 0.41 ±0.09 



2t 




FIG. 4: Plot of Tdecoh corresponding to upward transitions 
1 -> 2 vs R and P for 70 = 2.64. 

(weak coupling) and Tdecoh — 0.17 ± 0.02 (strong cou- 
pling). In order for the Markovian approximation to be 
valid, the decoherence time must be short compared to 
the characteristic decay times of the correlation function 
that determines the rate constant. 

Simulation of the master equation requires knowledge 
of the transition rates m a fs(Xo). These quantities were 
obtained by numerically integrating the time dependent 
memory function discussed above. In this calculation one 
must ensure that for a given Xo the transition is allowed. 
Otherwise m a p(Xo) is assigned a value of zero for that 
choice of subsystem coordinates. These restrictions were 
discussed in Sec. Mil This process is repeated for a range 
of Xo values generating the surface, m a p{Ro, Pq). We 
obtain a different surface for each transition (see Fig. O. 

The structure of these transition-rate surfaces is due 
entirely to classical evolution of Xo along the mean sur- 
face. It is precisely this evolution that leads to spread in 
the ensemble of trajectories giving rise to decoherence. 
Thus, even though the evolution we are ultimately in- 
terested in calculating is entirely in diagonal space, the 
probability of the nonadiabatic transitions is calculated 
from the off-diagonal or coherent evolution segments de- 
pendent on Xq. In this way decoherence is accounted for 
in the formalism. 



Simulation of the master equation 

Once the surfaces, m a p(Xo) are obtained, Eq. (pTTj) is 
simulated using the sequential short-time propagation al- 
gorithm where the probabilities of nonadiabatic transi- 
tions are given by II = \m a p\At/(l + |m a ^|Ai)4£ Note 
that the value of II is determined at each time step using 
the value of m a ^(Xo) corresponding to the specific value 
of Xo at that time. The initial sampling is taken from 
the spectral density function where the bath distribution 
is given by the conditional density, pjj(Xf,; Ro). 
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b) 




FIG. 5: Plots of m a p{Ro,Po) versus Ro and Po for 70 



2.64. (a) m\2(Ro, Po) portions of the surface have a value of 
zero, corresponding to regions where transitions are forbid- 
den due to insufficient kinetic energy in the subsystem. |(b)| 
rri2i(Ro, Po) does not have this feature as it corresponds to 
downward transitions where the subsystem gains kinetic en- 
ergy. 



The results of the calculation of the forward time de- 
pendent rate coefficient, kAB{t), are shown in Fig. [SI The 
figure compares the rate coefficients using adiabatic, mas- 
ter equation, and quantum-classical Liouvillc dynamics. 
As expected the plots show rapid decay on a time scale 
T m ic to a plateau region characterized by a much slower 
decay on the macroscopic chemical relaxation time scale, 
Tchem ~ 67 for weak coupling and w 1.3 x 10 6 for strong 
coupling^ 

In Fig. [HI we see that the short-time decay portion of 
the rate coefficient given by the master equation simu- 
lation is in agreement with the quantum-classical result. 
The time scale of this decay, T TO j C w 4 in the weak cou- 
pling case and rj 2.5 for strong coupling, is about one 
order of magnitude larger than the average decoherence 
time Tdecoh ~ 0.41 for weak coupling, and w 0.17 for 
strong coupling as discussed above. From the figures we 
conclude that indeed T de coh < r mic < T chem . This in- 



equality provides the conditions for the applicability of 
the Markovian approximation used to derive the master 
equation. The plateau regions for both quantum-classical 
Liouvile and master equation dynamics have lower values 
than those for adiabatic dynamics. The smaller rate con- 
stant for nonadiabatic dynamics is due to enhanced bar- 
rier recrossing as a result of motion on either the excited 
state or mean surfaces. The plateau value using master 
equation dynamics in the strong coupling case is slightly 
higher than that obtained using quantum-classical Liou- 
ville dynamics. This likely arises from the fact that in 
quantum-classical Liouville dynamics the system evolves 
on the mean surface for long times, allowing trajecto- 
ries to re-enter the region of high nonadiabatic coupling 
where quantum transitions take place. Thus, the rate 
coefficient is reduced due to recrossings in the barrier re- 
gion. In general, for both weak and strong coupling, the 
master equation provides quite a good description of the 
rate coefficient data. 



V. CONCLUSION 

The master equation calculations presented above bear 
many similarities to surface-hopping schemes that have 
been used previously to simulate nonadiabatic dynam- 
ics of quantum-classical system s 49 i 50 ' 51 i 52 ' 53 i 54 . It is use- 
ful to comment on some of the similarities and highlight 
the differences. In most surface-hopping schemes, and in 
our master equation dynamics, the classical degrees of 
freedom evolve on single adiabatic potential energy sur- 
face segments according to Newton's equations of motion 
governed by Hcllmann-Feynman forces. This dynamics 
should be contrasted with the trajectory evolution in 
quantum-classical Liouville dynamics, where the trajec- 
tory segments of the classical degrees of freedom evolve 
on single adiabatic surfaces as well as mean surfaces. 
The differences between our master equation dynamics 
and other surface-hopping methods lie in the prescrip- 
tion for quantum transitions and the manner in which 
decoherence is incorporated into the theory. For exam- 
ple, in the fewest-switches surface hopping scheme ; 49 ' 50 ! 51 
the probability of a transition depends on the nonadia- 
batic coupling matrix elements and the off-diagonal el- 
ements of the density matrix. In our master equation 
the probabilities of quantum transitions are determined 
by a Monte Carlo sampling based on the magnitudes of 
our phase space dependent transition rates, m a p(Xo), 
and the sampling algorithm reweights averages so that 
no bias is introduced. Decoherence is accounted for in 
the fewest-switches simulation by collapsing the density 
matrix onto a diagonal state depending on certain con- 
ditions such as motion outside a window of strong cou- 
pling^ In our calculation decoherence effects have been 
incorporated into the calculation of the transition rates. 

Decoherence has also been incorporated into the for- 
mulations of surface-hopping methods using other phys- 
ical principles. The idea in such methods is to include 
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a) 





quantum-classical Liouville dynamics^ and for the tem- 
peratures used in our calculations, these effects are very 
small. Regardless of whether the bath is treated clas- 
sically or quantum mechanically, decoherence enters our 
master equation through the forms of the transition rates 
and not as an additional term in the equation of motion. 

Finally, we remark that the simulation scheme for mas- 
ter equation dynamics has a number of attractive fea- 
tures when compared to quantum-classical Liouville dy- 
namics. The solution of the master equation consists of 
two numerically simple parts. The first is the compu- 
tation of the memory function which involves adiabatic 
evolution along mean surfaces. Once the transition rates 
are known as a function of the subsystem coordinates, 
the sequential short-time propagation algorithm may be 
used to evolve the observable or density. Since the dy- 
namics is restricted to single adiabatic surfaces, no phase 
factors enter the calculation increasing the stability of 
the algorithm. For complex reaction coordinates which 
are arbitrary functions of the bath coordinates the cal- 
culation of the transition rates will be more difficult and 
time consuming. Future research will determine if the 
master equation can be applied easily to realistic general 
many-body systems. Nevertheless, the results reported 
in this paper have served to provide a basis for an un- 
derstanding of the domain of validity of master equation 
approaches to quantum-classical nonadiabatic dynamics 
based on decoherence. 
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APPENDIX A: 



REDUCTION TO MEMORY 
FUNCTION 



the effects of decoherence in single trajectories, much like 
the description our master equation provides. For exam- 
ple, in the methods developed by Rossky et. a h 53 ' 54 i 55 
decoherence is introduced through an additional term 
in the evolution equation that accounts for the quan- 
tum dispersion about each classical phase space coordi- 
nate in the bath. In this way each of the trajectories 
in the full phase space experiences decoherence. The 
quantum dispersion of the bath is not included in the 
model calculations presented above but it is easily ac- 
counted for in our theory. In our master equation de- 
coherence arises through averaging over the bath phase 
space variables which were taken to be classically dis- 
tributed for high temperatures. The rate coefficient for- 
malism in Eq. (|29[) involves sampling from the full quan- 
tum spectral density function thus incorporating quan- 
tum dispersion in the bath coordinates. Such quantum 
effects have already been investigated in the context of 



Starting from the form of the memory kernel opera- 
tor given in Eq. (|10p , we may reduce this operator to a 
function. Without loss of generality, we assume that the 
adiabatic basis is real so that, C a p = C*^. Furthermore, 
taking the definition of J ^ and acting with the operator 
jav(Xo) coming from the leftmost J d '° operator on all 
operators to its right, the memory kernel operator may 
be written as 

M a p{t) 

= J2 D av (X )2Ke \u° av3v , (X, t) + U° avy0 {X, t) 
*D„,p{X Q )j v , p {Xv) ]av {X Q ) . (Al) 

In the above expression we have introduced the off- 
diagonal propagator, U° u ^ v ,{X,t) = (e" ir W) av ^ v ,. 
In the case of a two level system this propagator is given 
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exactly by 



(A2) 



for p^p'.u, v 1 = 1,2 and /i ^ /A ^ 7^ Here we used 
the fact that22i 

e (-^f' _i ' L M(*')( t_i ') = e - *-'* dT "W^ommV^- 1 ^' (*"*') 
= W wl /(t,*')e- <£ MM'(*-*') , (A3) 

to express the operator as a product of a phase factor 
and a classical propagator. We note that the only off- 
diagonal propagator matrix elements that contribute to 
the dynamics here are U° 2 12 — M$* 21- 

Recall from the definition of the momentum shift oper- 
ator, Eq. (fTij]) . that transitions can only occur if there is 
sufficient momentum in the subsystem to make a transi- 
tion to or from a mean surface. Otherwise the transitions 
are not allowed. Using the above form of the off-diagonal 
propagator in Eq. (|A1|) . the action of the memory kernel 
operator on some arbitrary function of the phase space 
variables, f(Xo,X b ), takes the following form: 

M al3 (t)f(X ,X b ) = 

S af3 2 2Re \w av (t, 0)1 D au {X 0nai )e- a "»'^ t 



x D l , a (XQ av )j va (X 0aiy )f(XQ aU7 X b ) 



+2Rc 



W a0 (t,O) D af s(X 0al3 )e 



(A4) 



xD a p(X 0a p)j a p(X 0a i3)f(X 0a i3, Xb) 



The arguments of / reflect the fact that we have acted 
with the rightmost momentum shift operator on the func- 
tion. If we now consider the action of the classical prop- 
agators that enter in the memory kernel operator we ob- 
tain, 

e iLa "^ Xa "^ t D ua (Xo au )j va (Xoav)f{Xo a i;, Xb) 

= D va (Xouv,t)jua(Xo a v t t)f(Xo au ji Xb : t) 

= D va {X 0au , t )f{XZ, t ,X b ,t) ■ (A5) 

In the last line we denoted the indices coming from the 
action of the second momentum shift operator as super- 
scripts {j va (X 0av )f(X 0av ,X b ) = f(X^,Xb)). Substi- 
tuting Eq. (|A5|) in the expression (|A4[) for the memory 
kernel we obtain, 

M aP {X, t)f(X , X b ) = 5 aP V M£(X, t)f(X™, X b , t ) 



the notation in the following calculation, it is convenient 
to write the generalized master equation (116[) in a more 
formal and compact form. Letting, 

j2M:$(x,t') P °i(x^ tn x b , t ,,t-t') 



+ £ M Z (x, t') P a d {xz v , v , x b ,t> ,t-t') 

V 

= (M(X,t')p d (X t ,,t-t')) a . (Bl) 
we can write Eq. p^|) as 
d 



Pd(X,t) = -iL d p d (X,t) 

+ f dt'M(X,t')p d (X t ,t-t') . (B2) 
Jo 



Of. 



Starting from Eq. (|B2p . we use standard projection 
operator methods 39 to obtain the evolution equation for 
the subsystem density matrix. If we let pg(A{,;i? ) be 
the bath equilibrium density matrix conditional on the 
configuration of the directly coupled i?o subsystem co- 
ordinates, we may define the projection operator as in 
Eq. ((201) and it's complement, by Q = 1 - V. Note that 
Vp d {X,t) = p c b (X b ;Ro)ps(Xo,t). 

Applying these projectors to the generalized master 
equation (|B2|) we obtain, 



0_ 

Of 



Vp d (X,t) = -ViL d Vp d {X,t) - ViL d Q Pd {X,t) 



+ / dfVM(X,t')Vp d (X t ,,t-t') 



dt'VM(X,t')Qp d (X t ,,t-t') , 



(B3) 



Of 



Q Pd (X,t) = ~QiL d P Pd (X,t) - QiL d Q Pd (X,t) 



+ / dt'QM(X,t')Vp d {X t <,t-t') 



dt'QM(X,t')Q Pd (X t ,,t-t') 



(B4) 



+Mf p (X,t)f(Xf 0v X b , t ), 
where the definition of M is given in Eq. (fTT)) . 



(A6) 



APPENDIX B: SUBSYSTEM MASTER 
EQUATION 

In this appendix we focus on the equation of motion 
for the subsystem density matrix. In order to simplify 



Solving the second equation formally we obtain, 

Q Pd (X,t) = e- lQL ^Qp d (X,0) (B5) 

- f dt'e- iQL ^'iQL d V Pd {X, t-t') + t) , 
Jo 

where the function 4>(X, t) involves fluctuations of the 
memory function from its bath average, 5M = M —{M)b, 
at various time displaced coordinates. Substituting this 
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solution into Eq. (|B3|) gives, 
d 

—p s (X ,t) = 

-(iL d ) hPs (X ,t) - J dX b iL d e- lQL ^Q Pd (X,0) 

+ f dt'(tL d e- lQL ^iQL d ) bPs (X,t-t') 
Jo 

+ f dt'{M(X,t')) b p s (X t ,,t-t') (B6) 
Jo 

+ ff dX b dM(X,t')e- iQL ^-^Qp d (X t ,,0) 

+ / / dt'dt"(SM(X,t')e- %QLdt iQL d ) b 
Jo Jo 

xp s (X t „,t-t' -t") 

+ Jo dt ' I dXbSM ^ X ' t '^ Xt '' t ~ t '^ ■ 

The last three terms in this equation involve integrals 
over the bath of expressions containing fluctuations of 
the memory kernel from its bath average. These expres- 



sions consist of SM correlated with dynamical quantities 
evolved under projected dynamics. By definition, SM is 
initially zero and, due to the presence of the phase factor, 
it oscillates strongly for long times. Consequently, the 
bath integral of the product of the oscillatory function 
SM with a time evolved dynamical quantity is expected 
to be small. Taking these considerations into account, 
we neglect the last three terms in Eq. (|B6[) . Making this 
approximation, the subsystem evolution equation takes 
the form, 

J^«tM)= (B7) 

- j dX b iL a e- lQLdt Q Pd (X,0) - (iL d ) bPs (X,t) 

- f dt'i(L d e- lQL ^'iQL d ) b p s (X,t-t') 
Jo 

+ f dt'(M(X,t')) b p s (X t ,,t-t') . 
Jo 

This equation, written explicitly in terms of its compo- 
nents is given in Eq. (f2Tj) and forms the basis for the 
reduction to a master equation. 
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